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The day-to-day 
changes in power 
and mechanical 
loads gets 
controllers out of 
tune. Ken shows us 
the math and 
algorithms behind a 
self-tuning controller 
that knows how to 
maintain efficiency. 
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^ very control 
HT engineer is familiar 

with the problem of 
tuning the parameters of 
a digital controller to suit the charac- 
teristics of a specific plant. Parameters 
perfect in one system may be wildly 
unstable or ineffective in another. 

The engineer adjusts the control 
parameters by trial and error — guided 
by instinct and experience, sometimes 
using heuristic rules — until the plant 
exhibits the desired behavior. The 
controller is then sent to the plant 
with the assumption that it will re- 
main reasonably constant. 

But, what if that assumption turns 
out to be false? Many systems have 
characteristics that change with power 
supply, mechanical load, or wear on 
the parts. In these cases, the quality of 
the control may degrade because the 
controller is no longer properly tuned. 

We need a self-tuning controller that 
automatically deduces a plant's charac- 
teristics and adjusts control parameters 
to maintain the desired behavior. 

The self-tuning control algorithm I 
present starts with a discrete-time 
proportional-derivative (PD) control 
algorithm. Piggybacked on that is a 
recursive least-squares (RLS) algorithm 
borrowed from adaptive DSP theory, 
which estimates the plant's character- 
istics from its I/O datastreams. 



The PD and RLS algorithms are 
well known and understood. What's 
not obvious is how to feed the results 
of the RLS back into the PD. 

SYSTEM ANALYSIS 

Figure 1 shows a typical system. 
Signal d{k) is the desired plant output, 
and y[k) is the actual output. The erroi 
signal x[k) is the difference between 
these outputs. The error signal is usu- 
ally represented as e(k), but I'm reserv- 
ing that symbol for something else. 
The controller output is u[k). 

For all signals, I prefer to normalize 
full scale to 1 . That way, the control 
algorithm can run independent of the 
application. Other software in the 
controller can include scaling factors 
as necessary. 

In discrete time, the plant output is 
described by: 

y(k + l)=-ay(k) + bu(k) (1 

where a and b are the plant param- 
eters. The controller function is de- 
scribed by: 

x(k)-d(k)-y(k) < 2 
u(k + l) = u(k) + K p x(k) + K d [x(k)-x(k-l) 

where K p and K d are the respective 
proportional and derivative control 
parameters. 

Assuming the controller can main- 
tain this plant at all, the output y[k) 
eventually reaches a steady state in 
response to a step input in d[k), at 
which time xfk) = within some band 
of tolerance. But, "eventually" doesn't 
necessarily constitute good, predict- 
able control. 

I'd like to optimize the plant's tran- 
sient response to a step input accord- 
ing to a defined measure by specifying 
the desired response time and over- 
shoot. 

In general, analytical solutions to 
discrete-time equations aren't possible 
so to describe the transient response, 



Plant k - 



Figure 1— In a typical feedback control system, the 
target and actual plant outputs are d(k) and y(k), 
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Listing 1— This listing shows the combined RLS-PD control algorithm in C. 



. include <stdio.h> 
^include <stdlib.h> 

/* Ruri-time Signals */ 
float "y, yzl; 
float u, uzl; 
float d; 
float x, xzl; 

/* Controller */ 
float Kp, Kd; 
float M, coeff; 



/* plant output y(k), y(k - 1) */ 
/* control signal u(k), u(k - 1) */ 
/* desired (target) plant output */ 
f* rate error x( k) , x( k- 1 ) */ 



/* PO control parameters */ 

/* intermediate value for Kp and Kd */ 



1 



, p21, 
1000 



p22 



/* RLS */ 
float 11. 12 
float gl. g2; 
float pi 1 . pl2 
#define I NIT_P 
float a, b; 
float alpha; 
int rls_cntr; 
#define RLS_CNT_VAL 10 
float get_actual_rate(void) ; 
float get_desi red_rate(void) 
void set_output(float) ; 
void rl s_pd_i ni t ( voi d ) 
( 

uzl - yzl = xzl = 0; /* 
rls_cntr - RLS_CNT_VAL ; /* 
a = b = 0; 
coeff =64/49 
Kp = 1: 
Kd = 1; 



intermediate values for vector g */ 

elements of vector g */ 

elements of matrix P */ 

estimate plant parameters, and */ 

elements of vector w */ 

estimation error */ 

RLS iteration counter */ 

/* called function prototypes */ 



clear run-time signals */ 
initialize RLS */ 



/* set arbitrary PD parameters */ 



.oid rls_pd(void) 

( /* Eq. 2, Control Output 

y = get_actual_rate( ) ; /* get run-time data */ 
d = get_desi red_rate( ) ; 



/ 



/* calculate rate error signal 
xzl); /* calculate control signal */ 
/* bound u to to 100% */ 



x = d - y: 
u += Kp*x + Kd*(x 
if (u > 1.00) 

u - 1.00; 
else if (u < 0) 

u = 0; 

set_output(u) ; /* set control output */ 

if (rls_cntr == RLS_CNT_VAL) { /* if time to reinitialize RLS */ 

pll = p22 = INIT_P; 

pl2 = p21 = 0; 



if ((--rls_cntr) == 0) { /< 
rls_cntr = RLS_CNT_VAL; 
M = (a + 1) / b; 
Kp = (64/49)*M*(a+l) ; 
Kd = M / 7; 



if time to calculate Kp and Kd */ 



11 - pl2*uzl - 

12 = p22*uzl - 
1=1+ i2*uzl 
gl = il / 1 : 
g2 = i2 / 1 ; 
alpha = y + a*yzl 
a += al pha*gl ; 

b += alpha*g2; 
pll -= gl*il; 
pl2 -= g 1 * i 2 
p21 -= g2*il 
p22 -= g2*i2 
xzl = x ; 
yzl - y; 
uzl = u ; 



pll*yzl;/ v 
p21*yzl ; 
- 11*yzl; 



Eq. 19, Gain vector g */ 



b*uzl; 
/* 
/* 
/* 



/* Eq. 22, Estimation Error alpha 

Eq. 21, Estimated Plant Parameters 

and b */ 

Eq. 20, Inv Correlation Matrix P * 



/* z A (-l), Increment Time */ 



I'll convert equations (1) and (2) to 
continuous time using these approxi- 
mations: 

f d (k) = f c (kx) (3) 

f c (t)4[f d (t + T)-f d (t)] 

= I[f d (t)-f d (t-T)] (4) 

where is the discrete-time-sampled 
version of the continuous-time vari- 
able /<., and x is the sample period. 

This approximation is not without 
risk. If the sample period is too long 
relative to the system's time constants 
or rise time, the first-derivative ap- 
proximation is too inaccurate to be 
useful. 

I developed the algorithm on a sys- 
tem with time constants 10-20 times 
the sample period. In this case, the 
magnitude and phase distortion intro- 
duced by sampling was negligible. 

The next step is to convert the 
continuous-time versions of equations 
(1) and (2) to a single equation describ- 
ing the response to a step input. By 
subtracting y{k) and dividing by x, 
equation (1) can be rearranged as: 

i[y(k + l)-y(k)] = 

-(a+i) yM + b tt(k) (5) 

Using the continuous-time approxi- 
mations in equations (3) and (4), equa- 
tion (5) becomes: 



y(t)=-l 



a+ 1 



'y(t) + fu(t) (6) 



If we identify: 

I 
J 

1 
J 



b 

' x 

a+ 1 



(7) 



then equation (6) simplifies to: 

Jy(t)-u(t)-Fy(t) (8) 

which is the familiar equation for a 
rotating system relating the moment 
of inertia (/), the coefficient of friction 
(F), and the driving force (u) to rota- 
tional speed (y). 

To derive the continuous-time 
approximation to equation (2), first 
define the continuous-time propor- 
tional control parameter K pc = K p /x. 
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Then, by rearranging and dividing by t, 
equation (2) becomes: 

I[u(k + l)-u(k)j- 

K pc x(k) + ^-[x(k)-x(k-l)] I 9 ' 

Using equations (3) and (4), equa- 
tion (9) becomes: 

u(t)-K pc x(t) + K d x(t) or 

u(t)-K pc d(t)-K 1> , : y(t) + K,,d(t)-K, 1 y(t) (10) 

By differentiating equation (8) and 
substituting (10), we obtain the sec- 
ond-order differential equation for the 
combined plant-controller system: 

Jy(t) + (K d + F)y(t) + K pc y(t) = 
K d d(t) + K pc d(t) 

The transient response of equation 
(11) to a step input needs to be ana- 
lyzed to quantify the response time 
and maximum overshoot. To analyze 
the transient response, take d(t) to be a 
unit step function and assume y(0) = 0, 
y'(0) = 0, and d(0) = 0. 

Then, taking the LaPlace transform 
of equation (11) and substituting: 



D(s)=4 

which is the LaPlace transform for a 
step function of 1, we get: 

Y(s ] K d S + K pc 

* W s(JsMK a+ F)s + K pc ) (12) 

The interesting characteristics of 
equation (12) can be gleaned by rewrit- 
ing it as: 

where the natural frequency (oj n ), the 
zero (-c), and the damping coefficient 
(£) are defined as: 

K d (14) 
K d+ F 

For the controller to be self-tuning, 
it must know what constitutes being 
in tune. 



The desired behavior of the plant 
must be defined in terms of a set of 
mathematical criteria that use data 
available to the controller. 

Three criteria to consider are the 
type of response (underdamped, criti- 
cally damped, or overdamped), the 
damping coefficient, and the maxi- 
mum overshoot to a step input. 

If some overshoot can be tolerated, 
the quickest convergence of the re- 
sponse to the desired value can be 
achieved with an underdamped sys- 
tem, meaning the damping coefficient 
is between and 1 . 

A damping coefficient close to 
yields a shorter rise time but greater 
overshoot. With a damping coefficient 
close to 1, overshoot is small but the 
rise time is longer. 

Unless there are other system re- 
quirements that affect the selection of 
the damping coefficient, 0.5 is a good 
compromise between reducing over- 
shoot and shortening the response. 

For a given damping coefficient in 
an underdamped system, one more 
factor affects the overshoot: the ratio 
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between the zero (-c) and the real part 
(-fa) n ) of the pair of complex poles [1]. 
If the ratio: 



is close to 1, the overshoot can be as 
high as 70%. Larger ratios reduce the 
overshoot, and by the time a = 4, the 
overshoot is down to about 20%. 

The minimum possible overshoot is 
about 14% as a approaches infinite. I 
used a = 16 to keep overshoot below 
20% allowing for some noise, but any 
number greater than 4 is just as good. 

By setting the damping coefficient 
(£) to 0.5, a = 16 reduces to: 



s/ K P c J =8K d 

From this and the definition of f in 
equation (14), it follows that: 



K d =, 
K J 

K p = T 



64 F z 



;49j) 

64F 2 
(49J) 



(15) 



RLS ALGORITHM 

Now that we can calculate K p and 
K d from the plant characteristics / and 
F, what next? We know the sample 
period x, and / and F are functions of a 
and b, in equation (6). 

So, if we can calculate a and b, we 
are done. To do that, I returned to the 
domain of discrete-time and borrowed 
the RLS algorithm from DSP theory. 

Least-squares is a method of finding 
the best solution to an overdetermined 
set of simultaneous equations. The 
solution to the least-squares equation 
is computationally expensive, involv- 
ing n-dimensional matrix inversion 
and multiplication, where n is the 
number of simultaneous equations. 

Such a calculation may be beyond 
the capability of a simple microproces- 
sor. But, there's a recursive solution to 
the equation in which a simpler set of 
calculations is performed for each 
equation individually, and the results 
from one set of calculations are fed 
into the next. 

To get the recursive solution, I need 
to cast the problem into a standard 



form known as the deterministic nor- 
mal equation. First define the data 
vector a: 

a T (k) = [-y(k-l),u(k-l)] (16) 

and the estimated parameter vector w. 

w T =[a,b] 

The performance measure used to 
judge the estimated parameters is the 
difference between the measured out- 
put y{k) and the predicted output ob- 
tained from the estimated parameter 
vector in equation (1). This estimation 
error e[k) is: 

e(k)-y(k)-(-ay(k-l ) + bu(k-l)) 
= y(k)-a T (k)w 

For a datastream of n samples, we 
get n instances of equation (16), which 
can be cast into matrix form by defin- 
ing the data matrix A: 



A(k. 



a T ( k - n + 1 ) 



(17) 
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Figure 2— Control is sluggish lor the 
first ten samples until new control 
parameters are calculated. At the 
second step input in d(kj, the 
controller is tuned. 



the error vector e: 

e T (k]=[e(k),...,e(k-n + l)] 

and the measured data vector y: 

y T (k)=[y(k),...,y(k-n+l)] 

From these definitions, e equates to: 

e(k) = y(k)-a(k)w 

In the least-squares method, the 
square of e should be minimized. The 
squared error function is: 

E(k)=e T (k)e(k) 

=y T y-y T Aw- w T A T y+w T A T Aw 

The surface defined by this function 
forms a bowl shape over the a-b plane. 
The a-b coordinates under the lowest 
point in the bowl define the point of 
minimum error between the measured 
output and the estimated output calcu- 
lated from the estimated a and b. 

At this point, the gradient of the 
error surface equals zero or: 

|§ =-2A T y + 2A T Aw 
= 

which reduces to: 



A T y = A T Aw 



(18) 



This is the deterministic normal 
equation for least-squares parameter 
estimation, and the solution w pro- 
vides the necessary a and b. 

Deriving the recursive solution to 
equation (18) is a rocky trek through 
linear algebra. But, now that the prob- 
lem has been made to fit Procrustes- 
wise into the standard form [2], there's 
no need to reinvent the wheel. Instead, 



I'll simply assume that a miracle oc- 
curs and jump right to the solution. 

First, I need to define one new ma- 
trix, two vectors, and a scalar. The 
gain vector g is defined as: 

K(k) _ P(k-Da(k) 

81 ' l+a T (k)P(k-l)a(k) lXy) 

where the inverse correlation matrix P 

is: 

P(k) = P(k-i)-g(k)a T (k)P(k-l ) (20) 

The estimated parameter vector w = 
[a,bj T is evaluated by: 

w(k) = w(k-l ) + g(k)a(k) (21) 

where the estimation error a is: 

a(k)=y(k)-a T (k)w(k-l) 

= y(k)-w T (k-l )a(k) (22) 

The RLS parameter estimation 
algorithm proceeds as follows. First, 
equation (19) calculates the gain vector 
g, and equation (22) calculates the 
estimation error (a). 

Using g and a, the estimated pa- 
rameter vector w is recursively up- 
dated in equation (21). In the last step, 
equation (20) recursively updates P. 
The algorithm is initialized at k = 
with w(0) = 0, and P(0) = pi, where p is 
an arbitrary number greater than 1. 



Once a and b are estimated, the 
control parameters K p and K d are four 
from equations (7) and (15), which 
reduce to: 



Figure 3— At first, the output wildly 
overshoots. Again, the controller is 
tuned at the second step input. 



K rf = 



7b 



(2. 



K f 64 ) (a + i; 
K p = l 49 J b~ 



In equation (23), the control param 
eters seem independent of the sample 
rate, and in a sense they are. The sam 
pie rate doesn't directly contribute to 
the control parameters, but a different 
sample rate would be made manifest 
by different plant parameters a and b. 

PD_RLS CONTROL ALGORITHM 

Several issues need to be considers 
before the RLS algorithm is ready. On< 
is the value of p used to initialize p(0) 

The only constraint on the value ol 
p — that the correlation matrix P~'(k) 
has to stay invertible — is met by set- 
ting p to a number greater than 1 [3]. 
Otherwise, p's value has little effect o 
the results, except that a larger num- 
ber results in faster convergence. 

Another consideration is when to 
recalculate the control parameters K p 
and K d . The estimated plant paramete 
vector w is inaccurate for a short time 
after initialization. Control parameter: 
calculated from immature parameter 
estimates don't result in good control. 

In a noisy system with a constant 
target rate, the system reaches a stead' 
state with the error signal equal to 
zero plus a noise component. In this 
case, the RLS algorithm can get con- 
fused and calculate plant parameters 
that reflect the noise more than the 
true plant characteristics. 

You could skip recalculating the 
estimated plant parameters unless the 
target rate changes by an amount sig- 
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nificantly greater than the magnitude 
of the noise. 

I found empirically that I can obtain 
jOod results by running ten iterations 
of the RLS algorithm before using the 
estimated plant parameters to recalcu- 
late the control parameters. Then, I 
reinitialize the RLS algorithm and 
recalculate the control parameters 
after every ten iterations to keep up 
with changing plant characteristics. 

Listing 1 shows the combined PD- 
RLS control algorithm in C. The inter- 
mediate vector i used in calculating g 
is also used in the calculation of P. 
This works because f = PA in (19), so F 
= A T P, which is what (20) needs. 

SYSTEM SUCCESS 

Figures 2 and 3 show data obtained 
from trial runs on the same system. In 
each case, the target rate is 0.5 from 
time to time 2.5 s, after which it 
goes to 0. At 4 s, it returns to 0.5 again. 

In Figure 2, the control parameters 
start at K p = 0.001 and K d - 1. The 
result is a sluggish response for the 
first second until the control param- 
eters are recalculated. The response to 
che input pulse at 4 s reaches the tar- 
get rate within 1 s, overshoots by about 
10%, and then settles down. 

In Figure 3, the control parameters 
start at K p = 1 and K d = 0.001, which 
results in almost 50% overshoot in the 
first second. But again, the response to 
the input pulse at 4 s is well-behaved. 

COMPUTERS IN THE FIELD? 

This algorithm allows the control- 
ler to calculate automatically the pa- 
rameters needed to achieve the desired 
response. It produces a transient re- 
sponse with specific characteristics, 
but other characteristics can be cho- 
sen. Different combinations of the 
damping coefficient and the zero in 
equation (13) simply change the con- 
stants used in equation (20). 

I implemented this algorithm on a 
68HC11 controller with a 12-MHz 
clock. I was able to use a sample pe- 
riod of 100 ms with plenty of time 
available for other system functions. 

The controller was used on a fertil- 
izer spreader — literally, "in the field" — 
in which the fertilizer application rate 
varied with the speed of the truck 



through the field to apply a constant 
amount of fertilizer per unit area. 

To complicate the issue, the ma- 
chinery's response changed with pow- 
er-supply voltage, load, and usage as 
bearings and gears jammed with fertil- 
izer and dust. Also, the operators 
wanted to use the same controller on 
several different systems and to switch 
from one to another at any time with- 
out retuning the controller. 

The algorithm proved to be robust. 
The controller could connect to a new 
system and be tuned in 1-2 s. And, the 
tuning changed with plant characteris- 
tics, keeping response consistent. |j§ 

Ken Baker works as a senior design 
engineer at Guidant/CPI, a medical- 
device company. You may reach Ken 
at qc03660@stp03.guidant.com. 



SOFTWARE 



The algorithm in 68HC11 assembly 
code, available on the Circuit Cellar 
BBS, makes use of floating point 
functions available on the Internet 
at http://freeware.aus.sps.mot.com: 
80/freeweb/amcu_ndx.html# 
mcull. HC11FP11 Asm is freeware 
maintained by Motorola. 
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